Oscillatory shear flow in nanochannels via hybrid particle-continuum scheme 
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This paper provides an application of our hybrid continuum-particle scheme to liquids by con- 
sidering unsteady shear flows driven by wall oscillations in nano-slots. The particle region (P) 
adjacent to the wall, is described at the atomistic level by molecular dynamics, while the outer 
region (C), described by a continuum equation for the transversal momentum, is solved via a finite 
volume method. Both regions overlap in the "handshaking" region where a two-way coupling scheme 
(C— »P and P^C) is implemented. A protocol for the C— >P coupling was presented in a previous 
' '" paper [Delgado-Buscalioni and Coveney, Phys. Rev. E, 67, 046704 (2003)]; here we focus on the the 

P^C counterpart, that is, on the exchange of information from the microscopic to the continuum 
domain. We first show how to use the finite volume formalism to balance the momentum and mass 
fluxes across the P— >C surface in such a way that the continuity of velocity is ensured. Then we 
analyse the effect of the fluctuations of the (P) stress tensor on the stability and accuracy of the 
numerical scheme. This analysis yields a condition which establishes a restriction on the resolution 
of the flow field that can be described by the hybrid method, in terms of the maximum frequency 
and wavenumber of the flow. 
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I. INTRODUCTION 

> 

■ Many systems, including complex flows near interfaces, wetting, drop formation, melting, crystal growth from a 
I fluid phase, and moving interfaces of immiscible fluids or membranes, are governed by the dynamic interplay between 
t-H ■ the rapid atomistic processes occurring within a small localized region and the slow dynamics occurring within the 
bulk liquid. Such problems are too computationally expensive for any standard molecular dynamics simulation (MD). 
A promising alternative is to use hybrid particle/continuum algorithms which retain the atomistic description only 
where it is needed, and reduce the computational cost by solving the bulk flow by much faster continuum fluid 
dynamics methods. 

Hybrid particle/continuum algorithms for solids and gases were the first to be fully developed in the literature. 
As also happens in theoretical descriptions, a hybrid description of the liquid state is the most challenging one. The 
• general procedure is to connect the particle region (P), described by molecular dynamics, and the continuum region 
(C), described by continuum fluid dynamics (CFD), with a handshaking region comprised of two buffers: C— >P and 
P^C. Most of the difficulties encountered when developing a hybrid scheme for liquids arise at the C-^P buffer, 
' where the molecular dynamics has to be reconstructed so as to adhere to the prescriptions coming from the C-flow. 
Within the P— >C region the microscopic variables are coarse-grained to supply boundary conditions for the continuum 
domain, although within this domain the molecular dynamics are not directly altered and their motion is uniquely 
governed by the corresponding molecular interactions. 

The kind of information to be transferred at the handshaking region has been the subject of some discussion. The 
first hybrid schemes to appear in the literature [see Delgado-Buscalioni and Coveney and references therein] consid- 
ered the coupling of momentum in steady shear flows using schemes based on coupling-through- state; in particular on 
matching the velocity of P and C within the overlapping region. Flekk0y et al. ,4] introduced another coupling strat- 
egy based on the exchange of the fluxes of conserved quantities (from and to P and C) and they applied their scheme 
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to transversal momentum and mass transfer in steady (Poiseuille) flows. This work was subsequently generalized to 
include the energy exchange . More recently, Delgado-Buscalioni and Coveney introduced modifications into the 
flux-based C^P coupling to allow mass, momentum and energy transfers in time-dependent flows. In that work, the 
MD region is surrounded by a continuum description and the C— >P protocol gives rise to the generalized boundary 
conditions for the particle system which provide the correct (hydrodynamic) propagation and relaxation of shear, 
sound and heat waves across the MD region Q . The same study compares schemes based on coupling-through-fluxes 
and coupling-through-state in the simulation of a relaxing flow due to decaying pressure waves. The results showed 
that the imposition of the continuity of the thermodynamic variables alone does not provide the correct decay of the 
hydrodynamic modes and, for example, it yields a negative entropy production |3|]. On the contrary, when using the 
flux-scheme the correct thermohydrodynamic behaviour was obtained. 

We would like to stress that, although most of the applications that any hybrid method must address involve 
unsteady flows, the previously proposed hybrid schemes for liquids have barely even considered time-dependent sce- 
narios. One of the main purposes of this paper is to study the applicability of the flux-scheme in the simulation of 
time-dependent shear flows, ranging from moderate to high characteristic frequencies. To that end we consider fluid 
flow in a slot driven by the oscillatory motion of one of the walls in its own plane (Stokes problem). This set-up 
retains the essential features encountered in unsteady shear flows while being analytically solvable. Oscillatory shear 
flow is widely used in rhcological studies of complex fluids, such as polymer brushes (see Wijmans and Smit 6] for a 
recent review). These systems are good examples of typical applications of the hybrid scheme, which can treat the 
complex fluid region by MD and the momentum transfer from the bulk by continuum fluid dynamics (CFD). The 
hybrid set-up we are considering here can be also used in the study of flow induced by complex boundary conditions 
(such as rough surfaces) in nanotechnological processes, an active area of research with important applications in 
micro electromechanical systems (MEMS). For instance, Stroock et al. showed that bas-relief nano-structures 
can be introduced on the floor of a slot in order to induce the mixing of solutions in low Reynolds number flows in 
microchannels. 

Another objective of the present work is to study the effect of the microscopic fluctuations on the spatial and 
temporal resolution of the overall coupling scheme. To that end we quantify and analyze the averaging procedure 
needed in the particle measurements in order to achieve a given signal-to-noise ratio. Finally, we consider the velocity 
discontinuities that arise at the particle-continuum interface. We show that these can be controlled by using a 
hybrid velocity gradient in the imposition of the momentum flux on the continuum system at the particle-continuum 
interface. The idea of using hybrid gradients was first proposed by Flekk0y and Wagner in work addressing the hybrid 
description of the density field in Fickian diffusion Q . 

The rest of the paper is organized as follows. In Sec^we present the domain decomposition used in our model, give 
details on the molecular domain and provide definitions of the spatial and time averages used in this work. In Sec. IIIII 
we present the hybrid finite volume method used to inject the averaged microscopic flux into the continuum domain. 
The protocol for the C— >P coupling is explained in Sec. IIVI The hybrid scheme produces a finite signal-to-noise ratio 
due to the fluctuation of the particle momentum flux and, in Sec. E] we discuss how to use averaging procedures 
to maximize it. In Sec IVII we apply the hybrid scheme to compute oscillatory shear flow in nanochannels. We also 
present some hybrid simulations to verify the conclusions of Sec. [V] concerning the role of fluctuations. Concluding 
remarks are given in Sec. IVIII 

II. THE MOLECULAR REGION 

A typical spatial domain decomposition for our hybrid scheme is depicted in Fig. ^ In the present set-up the 
overall domain ranges from x = to x = L x and from — [L a /2, L a /2] in the other two periodic directions (a = {y, z}). 
The particle region (P) spans from x = to x = xqp, the continuum region (C) from x = xpe to x = L x and both 
subdomains overlap within xpe < x < xqp- 

The particle region (P) contains N(t) particles at time t, interacting through interparticle potentials and evolving 
in time through Newtonian dynamics. In this work, the particles interact through a truncated Lennard- Jones (LJ) 
potential ip(r) — ipL.j{r) — ipLj{r ) where tpLj(r) — 4e~ 1 [(a/r) 12 — (er/r) 6 ]. Simulations were performed for the 
purely repulsive WCA Lennard-Jones potential using a cut-off radius of r Q = 2 1 / 6 er and for the standard LJ fluid 
using a larger cut-off radius r Q = 3a. 

Each particle has a mass m, velocity Vj and energy = ^mvf + SjV-'( r y) ( r ij = r j ~ r i)- The equations of motion 
for the particles are 



i\; = Vj 



(1) 
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N(t) 

in = f i = ^2rr j 1 dip(n j )/dr ij T i:i , (2) 

i<j 

where fi is the the force acting on the i th particle. Equations Q and were solved via standard molecular 
dynamics (MD) using the velocity- Verlet algorithm with a time step Atp ~ 10~ 3 t, where t = (mcr 2 /e) 1 / 2 is the 
characteristic time of the LJ potential. Throughout the rest of the paper, all quantities will be expressed in reduced 
units of the LJ potential: r(= 0.45 x KT 13 s), a(= 3.305 x 10~ 12 cm), e, m(= 6.63 x 10 _23 g) and e/k B (= H9.18K) 
for time, length, energy, mass and temperature respectively (the numerical values correspond to argon). 

The particle system is bounded in the x direction by a lower wall located at x — which contains a layer of 
1600 spheres strongly tethered to the sites of the (1,1,1) plane of an fee lattice by harmonic springs of stiffness 
k = 1320eer~ 2 . The wall atoms do not interact with each other and the wall-fluid interaction is LJ with an increased 
cutoff and potential well, r o tu = 1.311cr and e w — 1.303. These values ensure vanishingly small slip velocities for the 
shear rates considered here. 

In order to maintain the temperature of the system constant we used either a Nose-Hoover thermostat ^l| for the 
LJ fluid, or a Langevin thermostat |TTL ll2T ] for the WCA-LJ fluid. The local kinetic temperature, needed for the 
Nose-Hoover thermostat, was calculated from the peculiar velocities within slices of width l.Ocr along the x direction. 
Also, to ensure that the Langevin thermostat does not bias the shear profile, the Gaussian white noise and damping 
terms are only added to the equations of motion for the velocity components normal to the mean flow, the x and z 
directions [l2| . Both thermostats proved to be equally efficient. 

The coupling scheme is a two-fold communication between C and P, as can be understood from Fig. ^ Within 
the C— >P domain, the fluxes evaluated from the continuum solution are imposed on the particle dynamics. This 
part of the coupling scheme has been discussed by Delgado-Buscalioni and Coveney for the general scenario of 
mass, momentum and energy transfer in unsteady relaxing flows. On the other hand, the flux due to the particle 
dynamics across the P— >C interface is to be injected into the continuum system as flux boundary conditions. In 
order to coherently communicate information between both descriptions of matter, the fluxes arising from the particle 
system must be averaged in time and space prior to transferring to the coarse-grained level. Moreover these averages 
need to be local on the coarse-grained coordinates and on the time scale of the C domain. To that end, around each 
coarse-grained spatial coordinate, R, we define a cell of volume ^(R) and, for any particle variable, say </>(rj,i), we 
define the local mean <&(!?, t) at the coarse grained location R as 



jV(R) £i=i <f>(ri,t)8(r - r,)d 3 r 



7 ; w - — - , (3) 



where the integration is performed over the cell volume V^(R) and the summation runs over the particle index. The 
corresponding spatio-temporal average is also local in the "coarse-grained" time tc with a time averaging window of 
length At av , 



($}(R,ic) - -rj- f ° $(R,t)6(t-t k )dt. (4) 

AX a v Jtc-Atav 



For the temporal average a number of samples n s of the microscopic quantity are taken along the time interval At av 
at times given by f& — tc — At av + kSt s (k = {1, ...,n s }). To ensure statistically independent measurements the 
time interval between samples St s cannot be smaller than the decorrelation time of <3 ) , which can be calculated from 
f™(&(t))&{0))dt/{{&) 2 ), where $' = $ - ($) (see Sec. 

As shown in Sec. II I II the average defined by Eqs. and |@J is used within the P^C cell to evaluate the 
coarse-grained momentum flux from the particle region to the continuum domain. In this case the spatial average 
is computed within the volume Vpc = V(Rpc), where Rpc is the location of the P^C cell. This average shall be 
denoted (<&) pc(t) = (®}(R-PC,t)- We shall focus on the P^C coupling of the shear stress. At any given time instant, 
the time-averaged mean momentum flux tensor evaluated within the P— >C cell is given by 

Up) = ^^(v^mviVi-^fjCrvF^-npc, (5) 

where Npc is the number of particles inside the P— >C cell, ppc = Npc /Vpc, and npc is the surface vector of the 
cell pointing toward the neighbouring C-region (see Fig. H 3 ). 
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III. CONTINUUM DESCRIPTION AND P^C COUPLING 



Within the C region the relevant variables are the macroscopic local densities associated with the conserved quan- 
tities. If $ is any conserved quantity (per unit mass), its conservation law is given by 

£~V.J., (0, 

In this work we shall consider isothermal fluids and restrict ourselves to the momentum density pu(r,t). In what 
follows, the continuum velocity will be denoted by u while (v) stands for the mean particle velocity. The flux of 
momentum is given by J u = puu + II, where the pressure tensor II = P 1 + r includes the local hydrostatic pressure 
P and the viscous stress tensor. For an LJ fluid, the latter satisfies a Newtonian constitutive relation [T^L IT3 |. 
t = —r) (Vm + Vm t ) + (277/3 — £) V • u. where 77 is the shear viscosity and £ the bulk viscosity. In principle, the set 
of equations given by Eq.© can be solved by any standard continuum fluid dynamics (CFD) method. Nevertheless, 
as the hybrid scheme proposed here is based on the balance of fluxes it fits naturally with the finite volume method 
[l5j . which exactly balances the fluxes across the computational cells. 

In order to illustrate the P— >C coupling protocol we shall consider the flow of an incompressible and isothermal 
fluid with mean density p c and mean temperature T c . The fluid fills the space between two parallel walls at x = 
and x = L x . The wall at x — L x is moving with arbitrary velocity u wa ii(t) along the y direction. The flow is uniquely 
driven by the motion of this wall, meaning that the mean pressure is constant throughout the domain, P = P c {p c , Tc), 
and that there are no transfers of mean energy or mass along the x direction (perpendicular to the P— >C interface). 
As discussed by Vogel et al. , if the system is not confined in the flow direction (as happens to be the case in our 
system) this flow is laminar and the velocity of the continuum flow is given by u = u y (x,t)j. The stress tensor is 
J u = P c 1 — i]du y /dxj. In the case of the LJ fluid, the value of the dynamic viscosity 77 was obtained from Heves Il4l 
and for the WCA-LJ fluid 77 was measured from MD simulations using a standard non-equilibrium procedure |17| . 
i.e., from the ratio of the shear stress and shear rate in Couette flow. The equation of motion for the y- momentum is 
given by the unsteady diffusion equation 

dUy _ C^Uy 

dt dx 2 ' {> 

where v = 77/ p c is the kinematic viscosity. At the fixed wall, the boundary condition is u y (0, t) = and at the moving 
wall u y (L x ,t) = u wa u(t). 

For the finite volume discretization of Eq. (|7| in the geometry of Fig. ^ we divide the extent of the C region 
along the ^-coordinate (xpe < x < L x ) into M cells. The faces of the cells can be placed at arbitrary positions 
x{ with k = {0, ...,M}, and the centres of the cells at Xk = {x[ + x{ 1 )/2 with fc G {1, ...,M}. To derive a closed 
set of equations for the velocities u y (xk) at the fc G {1, ..,M} cell centres, one integrates Eq. {7J) over the volume 
of each cell. For a given cell, say k = H, the integration is done over x w < x < x e , where the cell faces of the H 
cell (x w and x e ) are illustrated in Fig. ^p. The resulting first order spatial derivatives at the cell faces x e and x w 
are discretized as, e.g., du y (x e )/dx = (u v (xe) — u v (xh))/(xe — xh); where the "east" cell E = H + 1 is placed at 
xe = xh+i and the "west" cell W at xw = %h—i- We use an explicit time integration scheme by approximating the 
time derivative at the centre of the cell by du v /dt ~ (u y (t + Atc) — u y(t))/ Ate, where At c is the time step. We 
refer to the excellent book of Patankar |15| for a comprehensive presentation of the finite volume discretization of 
Eq. Q and Navier-Stokes equations including a discussion on the different time integration schemes. The resulting 
discretized version of equation reads, 

u y (x H , t + At c ) = (1 - r e - r w ) u y (x Hl t) + r e u y {x E , t) + r w u y (x w , t), (8) 
with E = H+1,W = H — 1 and H = {1, .., M} . 

We have introduced r e and r w , defined as 77 = vAt/(Axf Axr) where / = {w, e}, Ax? is the the centre-to-centre 
distance and Axh the face-to-face distance, as shown in Fig. lb. Most of the simulations reported in this work 
were performed with a regular grid 31] (7' = r w = r e ), using Ax ~ 0.5. The time step Atc was chosen to fulfil the 
Courant condition r < 1/2, which ensures the stability of the explicit scheme in Eq. JSJ) 15] . Note that Eq. JHJl is 
not closed for the velocities at the H — {1, .., M} cell centres. For instance, the evaluation of the velocity at the cell 
H = M requires the velocity at the outer cell E = M+l, which is outside the computational domain of the C region. 
Similarly, for the velocity at H = 1 one requires the velocity at the outer cell W = 0. The velocities at these outer 
cells (also called ghost cells) have to be evaluated from the boundary conditions. For instance, the velocity of the 
ghost cell E = Al + 1 at Xm+i — L x + Axe/2 is calculated by linear extrapolation from the velocity at the H = M 
cell and the velocity at the wall: u v {xe, t) = u y (xn, t) + 2Ax e {u wa ii{t) — u v (xh ,t))/ Axh , where H — M. 
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In the same way, to close Eq. J5J for H — 1 we need to determine the velocity at the ghost cell W — at 
x = xpc — Axw/2- To that end we use the hybrid formulation and set the y-momentum flux across P-^C equal 
to the time-averaged momentum flux due to the particle system at the P— >C cell: {j p , xy )pc — ~ri{du y /dx) Xpc . The 
discretized version of the right hand side (RHS) of this momentum flux balance equation is expressed as —rj(u y (xH) — 
u y (xw)) I Ax w , with H = 1 and W — 0. Note that both the particle and the continuum domains overlap at the first 
cell H — 1, so the velocity u v (xh) could be chosen to be either the time averaged mean particle velocity at the H = 1 
cell, {v v )(xh) or the continuum velocity at the same location, u v {xh)- The implication of such a choice will become 
clear shortly. We have enabled a set of possible choices based on a linear combination of both velocities, 

u v {xh) = (1 — a)u y {xH) + a ( v y){ x ii), for H = 1 and a G [0, 1]. (9) 

The flux balance equation (j PtXy )pc = ~ r l(^ L y( x H) ~ u y (x-\y)) / Ax w provides the velocity at the ghost cell u y (xw), 
required to close Eq. © for H = 1, 

u y {x w ) = u y (x H ) + (j p , xy )p C Ax w h for W = and H = 1. (10) 

Insertion of Eq. (|10|) and Eq. into Eq. JSJ provides the time advancement algorithm for the velocity at the cell 
interfacing the particle system: 



u y (x H ,t+At c ) = (1 - r e ) u y (x H , t) + r e u y (x E ,t) + \^y)pcif)At 

pAx w 



(11) 



+ ar w ({Vy)(xH,t) — u y (xH,t)) , for H =1,E = 2. 



The effect of the choice of a in Eq. © becomes clear by inspection of Eq. (|lll) . For any a / the last term 
on the RHS of Eq. 111(1 acts as a relaxation term that ensures velocity continuity by gently driving the continuum 
velocity at the boundary cell x — x± to the corresponding particle average (v y )(xi). This term vanishes or becomes 
negligible once velocity continuity is established, u y (xi) ~ (v y ){x\). We emphasise that this sort of velocity coupling 
does not directly act on the particle dynamics, but only within the algorithm for the continuum dynamics. This fact 
is significant because it avoids any artificial alteration (such as a Maxwell daemon) to the microscopic dynamics of the 
particles near the P— >C interface. The idea of using the scheme in Eq. I|ll|) arose from the outcome of calculations 
performed at low shear rates for which du y /dx < 10~ 2 . At low shear rates the velocity fluctuations become significant 
enough to produce deviations of the mean particle and continuum velocities at the overlapping region. If the velocity 
coupling term is absent in Eq. I|ll|) (a — 0) the deviations of the P and C velocities produced by subsequent particle 
fluctuations can be maintained in time, although the slopes of the P and C velocities do converge to the same value 
(see Figs. EJIHland discussion below). The reason for this fact is clear: the imposition of a flux does not prescribe the 
value of the velocity, but only its gradient. As shown below, the algorithm in Eq. (|llfl solves the problem of velocity 
continuity even when a small velocity coupling is used a (~ 0.1). 

To demonstrate the effect of the parameter a in Eq. lllfl on velocity continuity we performed some hybrid sim- 
ulations of Couette flow using a Lennard- Jones fluid with density p = 0.8 and temperature and T = 1.0. Similar 
tests were also done in the simulations of oscillatory shear flow presented in Sec IVII These tests were performed 
with a = {0,0.1,0.2, 1}. Figure |21 shows the y-velocity profile of a Couette flow with du y /dx — 0.01724, obtained by 
averaging within slices of width 0.2<r and over the whole simulation time. Results obtained for a > do not present 
any velocity discontinuity and are in perfect agreement with the analytical profile (the local shear rate at the xcp 
interface deviates by less than 2% from u wa u/ L x ). On the contrary, when using a = in Eq. I|ll|) . the averaged 
velocity profile presents a significant discontinuity at the overlapping region which affects the overall shear rate of 
the flow. As shown in Fig. [2 the discontinuity tends to be reduced as the average time is increased, but it is still 
significant after an average time over 1500r. 

The effect of the parameter a in the continuity of velocity is more clearly shown in Fig. [3] where the continuum 
and mean particle velocities at the boundary cell H = 1 are plotted versus time. Using a = 1 the relative differences 
between the C and P velocities are less than 2%, while for a value of a as small as 0.2, these relative differences are 
still less than 5%. By contrast, when using a = 0, this relative difference increases to ~ 80%. As seen in Fig. |3 
for a = the C velocity and averaged P velocity at the cell H — 1 "oscillate" around the correct value ~ 0.28, and 
the period of these oscillations can be as large as ~ 600r. For smaller shear rates the situation worsens because this 
period may become much larger. In summary, using a = the numerical scheme of Eq. Ijlll) does not guarantee 
velocity continuity. Inspection of Fig. 3a clearly indicates that using a > enables us to deal with time-dependent 
flows, as shown in Sec. IVII 

The central idea of our hybrid scheme is the balance of fluxes and such a balance must to be respected by any choice 
of a. To investigate this we compared the time averaged flux injected from the particle system with the time averaged 
flux measured from the continuum field at the interfacing cell, H = 1. If the flux exchange is correctly balanced, both 
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Particle 


system 


Continuum system 


a 


(jp,xy)pC 




n xy (xp C ) S n 





0.0289 


0.09 


0.0294 0.01 


0.2 


0.0287 


0.08 


0.0294 0.01 


1 


0.0298 


0.08 


0.0297 0.01 



TABLE I: The fluxes across the P^C interface for different values of the a parameter. Comparison is made between the time 
averaged particle flux (jp,xy) an d the time averaged momentum flux in the continuum field, Yl xy (xpc) = —r}(du y /dx)pc- Time 
averages were performed over ~ lOOOr. The results were obtained from a Couette flow with imposed shear rate u wa u/L x = 
0.0172 (see Fig. |2Jl. The momentum flux measured from the C flow is calculated using r\ — 1.75 ±0.05 and a discretized version 
of the velocity gradient near the P^C interface, du y /dx — (u(x2) — u(xi)) / Ax (where xi = xpc + Ax/2, Xi — x\ ± iAx and 
Ax = 0.9 is the width of the spatial mesh). For all values of a the time averaged P and C-fluxes coincide well within their 
fluctuation margins and agree with the externally imposed stress, rju wa u/L x — 0.030 ±0.001. The amplitude of the fluctuations 
of both fluxes (their standard deviation) has been also indicated by Sj and Sn- All quantities are expressed in LJ reduced 
units mo- _1 T -2 . 

quantities should coincide. The outcome of these tests is shown in Table 1. For any a E [0, 1], both fluxes differ by 
less than 2% indicating that the velocity coupling term in Eq. l|llf) does not alter the balance of fluxes. 

We note that the velocity coupling term in Eq. (|llf) introduces an extra momentum flux into the C system whose 
magnitude is arj ((v y )(xH,t) — Uy(xn,t)) /Ax. This flux fluctuates on the microscopic time scale and its average is 
close to zero. This fact explains why the relaxing term does not affect the flux balance on the time scale of the 
coarse-grained dynamics, as can be seen in Table 1. We measured the relative size of the flux contributions in Eq. 
(|ll|l and concluded that the size of the relaxing term remains much smaller than the particle flux for shear rates larger 
than 7 > 10~ 3 . Below a certain shear rate, 7 ~ 10~ 3 , the fluctuations of the relaxing term are of the same order as the 
particle flux signal. Nevertheless, for such small shear rates the signal-to- noise ratio of the particle flux becomes larger 
than one (see Sec. 0) and we are obliged to perform temporal averages over long times to recover the correct velocity 
profile and flux balance. We performed simulations of a Couette flow at 7 = 10~ 3 using a = 0.5. The values of the 
fluxes, averaged over a time of 10 4 r, were (j p , X y)pc — (1-9 ± 2.1) x 10~ 3 while —rj{du y / dx) pc = (1-8 ± 0.3) x 10~ 3 
(here the error bars denote half the standard deviation). These values compare well with the externally imposed 
stress, rju W aiil L x — (1.75 ± 0.05) x 10~ 3 . In passing, we note that owing to the effect of the strong fluctuations, 
when simulating flows with 7 ~ 10 -3 using a = 0, the averaged P and C velocities in the handshaking region never 
matched; hence, although the fluxes (and therefore the velocity gradients) were perfectly balanced, they did not agree 
with the desired, externally imposed, value. 

One might ask what value of a ^ is best suited for a certain simulation. We give here indications as to how to 
answer this question. Note that the velocity coupling takes a characteristic relaxation time 0[Ax 2 /(va)\ to drive the 
continuum velocity towards the average particle velocity. As stated above this time can be very short, for instance 
in simulations with a = 1, Ax ~ 1 and v ~ 2.0, one finds a relaxation time ~ It which lies within the characteristic 
time of the microscopic fluctuating currents |3j . In flows where the amplitude of particle fluctuations is expected to 
be large (e.g. at very low shear rates or near the critical point) it may be desirable to increase this relaxation time 
by using a smaller value of a ~ 0.1. In this way one can prevent large jumps in the continuum velocity which may 
lead to numerical instabilities of the CFD algorithm for the continuum domain. 

IV. THE C^P COUPLING: CONTROL OF MASS FLUCTUATIONS 

For a detailed description of this part of the hybrid scheme we refer to Delgado-Buscalioni and Coveney [3 , where 
the C— >P protocol is applied to the general scenario of flows with time-dependent mass, transversal and momentum 
and energy transfers. For completeness we now briefly explain how the momentum flux from the C flow is injected 
into the particle system. We also present a way to control the fluctuations of mass across the C— >P interface so as to 
ensure the balance of mass flux. 



A. Imposition of the shear stress on the particle system 

This part of the coupling scheme occurs within the C^P cell shown in Fig. The stress induced by the C-flow 
in the P domain is given by the local momentum flux at the C— >P interface at x = xcp- The flow considered here 
only carries transversal momentum (i.e. along the y-direction) , so the momentum flux is given by 



n C p ■ n C p = -P<A ± r]{du y /dx) XCP i, 



(12) 
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where ncp = — i is the surface vector of the C^P interface. In order to introduce this stress into the particle dynamics 
we add an overall external force F ext = AUcp ■ ricp to those particles within the C^P cell. For any instant of time, 
t, this force is equally distributed among the Ncp{t) particles inside the C^P cell, so the external force per particle 
is F ex t /Npc = AIIcp • ricp/Ncp- Note that this external force has a component normal to the C^P interface, 
which provides the hydrostatic pressure, and a tangential component providing the shear stress. The particles are free 
to enter or leave the C^P region, so the number of particles within this region Ncp(t) and the value of the overall 
external force fluctuate in time. The mean value of Nqp is AAxcpPcp so the average "pressure force" per particle 
is P c / (AxcpA 2 pep), where Axcp — 2a is the extent of the C— >P cell along the x direction, A — L y x L z , the local 
density is pep = Ncp/Vcp, and the pressure P c (p c ,T c ) is given from the equation of state which in the case of the 
WCA-LJ fluid was provided by Hess et al. and for the standard LJ fluid by Johnson et al 0. This normal force 
prevents the escape of particles and maintains the correct value of the density across the inner part of the MD domain 
as seen in Fig. 0] The shear force is distributed over the particles in the same way as described above for the pressure 
force. In this case, the flux of 2/-momentum to be injected in the particle system is rj{du y /dx) XCP . 

B. Control of mass fluctuations across the C^P interface 

For the flows we are treating here, the mean fluxes of mass and energy across the C and P interfaces are zero 
but fluctuations in the particle system produce perturbative mass currents along the x direction which need to be 
controlled. For this purpose, we coupled the end region of the P domain with a control equation that governs the 
mass flux across the particle system frontier. This control equation is derived using the finite volume rationale: we 
integrate the mass conservation equation dp/ dt = — V(pu) in a control volume around the C— >P interface (between 
xcp — xcp — Axcp/2 and xq = xcp + Axcp/2) to obtain the particle flux towards the C domain. Zero mass flux 
is ensured by equating this rate to the rate of insertion in the particle system, Npc- This yields 

N C p = A ((pv x ) C p - p c u x (x )) = A(pv x ) C p 1 (13) 

where Ncp is the time derivative of the number of particles within the C— >P cell. The mean particle mass flux 
{pv x )cp is evaluated within the C^P cell and averaged over At av . On the right hand side of Eq. I|13(l we have 
used the fact that the normal velocity in the C domain is zero, u x (xq) = 0. Equation proved to be successful in 
ensuring mean zero mass flux and in eliminating any incoming longitudinal mass currents towards the C— >P interface. 

We also used a different strategy to fix the rate of particle insertion/removal, which provides further control on the 
particle density near the C-^P interface. In this second approach the particle insertion rate is chosen to make the 
particle number in the C^P cell pc p relax to a pre-set value po , 

Ncp = — ((p)cp ~ Po) , (14) 

where Vpc is the volume of the C^P cell and r r is a relaxation time. As before, (.} denotes the time average over 
At av . The amplitude of the density fluctuations within the C^P cell is controlled by the coefficient r r and can be 
chosen accor ding to the requirements of the problem. For instance, in hybrid simulations of tethered polymers under 
Couette flow |l0t one may need to smooth out the perturbative density waves originated by the tumbling motion of 
the polymer at the inner part of the MD domain. Otherwise these waves can bounce back at the C-^P interface and 
affect the dynamics of the polymer. To that end, the relaxation time r r need to be smaller than the time needed 
for a sound wave to cross the width of the C— >P buffer, Axcp/c s , where c s is the sound velocity (see e.g. Q). The 
procedure implemented according to Eq. ifTIfl is an artifact that ensures fluctuations carrying mass and longitudinal 
currents are filtered out from the simulation box. In problems where the pressure waves are an essential part of the 
study, the mass flux will be governed by the longitudinal momentum and energy equations, which will need to be 
added into the hybrid formalism. According to Eq. (|14|) , particles are extracted if Npc < and, as explained in Ref. 
0, the first particles to be extracted are those closest to the C— >P interface. If Npc > 0, new particles are inserted 
with a velocity extracted from a Maxwellian distribution with mean velocity and temperature given by the continuum 
values u(xcp) and T — T c . The insertion of particles in dense fluids is a far from trivial task but it has been solved 
by the USHER algorithm proposed by Delgado-Buscalioni and Coveney |20| . 

The value of po in Ea. (|14f> was set to a slightly smaller value than the mean density. The reason for this choice 
was first to alleviate the computational cost of insertion and second to reduce the amplitude of ripples of the 
density profile near the C— >P buffer. For a fluid with p c = 0.8, we compared the density profile arising from a full 
MD simulation with that resulting using po = 0.7 and po = 0.5 in Eq. 114(1. This comparison is shown in Fig. 0] 
For both values of po, the "hybrid" density profile presents some ripples near the C— >P cell which are damped after 
around 3a. Within the P^C cell the hybrid density profile perfectly matches the density within the bulk. 
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We also estimated the cost in CPU time of the usher algorithm for the particle insertion (see Ref. 01 f° r further 
details). As an example, in a fluid with p c — 0.8, using p = 0.7 in Eq. lfH|) . a total number of particles of ~ 10 4 
and an overall simulation CPU time of 20000 seconds, the CPU time spent within the particle insertion/extraction 
subroutine was about 0.003 times the CPU time used by the main force subroutine of the MD code used in this study 
(which uses the Verlet list for counting neighbouring particles 10] ). This excellent performance indicates that the 
usher scheme is well suited for the insertion protocol in the hybrid scheme when using Lennard- Jones particles as 
the fluid of interest. 



V. SHEAR STRESS FLUCTUATIONS AND ACCURACY LIMITS 



In our scheme, the fluctuating nature of the fluxes introduced into the C region from P— >C imposes a limitation on 
our ability to resolve the flow field, as also arises in experiments and full MD simulations. This limit is determined by 
the signal-to-noise ratio becoming smaller than one. We shall now consider the effect of fluctuations of the transversal 
momentum flux and provide an estimation of its signal-to-noise ratio in order to derive the resolution limit of the 
hybrid scheme. This task requires consideration of the autocorrelation of the transversal momentum flux or, in other 
words, of the xy component of the microscopic stress tensor defined in Eq. (J5J 

G(t) = ^(f p , xy (t)f p , xy (0)) s u (15) 

where {.} s t denotes the statistical ensemble average (assuming ergodicity, we measured it by performing the time 
average over a long enough simulation time) and jL xy = jp,xy — (j P ,xy)st- Equation 1(15(1 gives a direct relationship 
between the amplitude of the stress tensor fluctuations and the shear modulus G(0) 

{(j' P ,xy) 2 )s t = ^ G ^)- ( 16 ) 



In the case of a LJ fluid the shear modulus is known and is given by G(0) = 3P - ^pU - 2pT |2l Ha]- For other 
interatomic potentials the route described above would require the evaluation of G(0) from MD simulations since the 
algebraic relationship of G(0) with thermodynamic variables is not normally available. 

To evaluate the size of the stress tensor fluctuations we shall use an equivalent route involving the transport 
coefficient (shear viscosity) and the corresponding Green-Kubo relation 

i]= G{t)dt. (17) 
Jo 

Let us approximate G(t) w rjexp(—t/TQ)/TG, where tq is a decorrelation time which can be calculated from the 
outcome of a simulation by r G = G(t)dt/G(0). Using Eqs. (|T5)l - ljT5l one obtains 



{(f P , X y) 2 )st = (18) 
Vp C T G 



We note that lim TG ^o[ ex P(~^/ T G)/ T G] — 25 (t), where 5(t) is the Dirac delta function. Hence the Landau expression 
for the fluctuation of the stress tensor II is coherently recovered from 1(18(1 in the limit of the delta correlated process, 

T~G > 

((Ky) 2 )st = ^T-m- (19) 

" Vpc 

The Landau expression p3| comes from a Langevin dynamics approach and is valid for a coarse-grained process in 
which the time step (or time between measurements) is much larger than the decorrelation time tq. 

Table 2 shows the amplitude of the fluctuations of the microscopic shear stress ((j' p xy ) 2 ) for several calculations of 
a Couette flow at very small shear rate 7 ~ 0.002. We confirmed that for these values of 7 the variance of the shear 
stress and velocity are similar to their values at equilibrium. Results were obtained for both the LJ and the WCA 
fluid. The numerical results are in good agreement with the theoretical expression given in Eq. 1)18(1 , where the decay 
time was measured from simulations, tq ~ 0.06. In the case of the LJ fluid (for which the value of G(0) is available) 
the values of (U' PtXy ) 2 ) obtained from Eq. ((TBI are compared with Eq. ((T5|> in Table 2. 

The standard deviation of the microscopic momentum flux averaged over a time interval At av = n s St s is given 
by S 2 — ((j'p xy ) 2 ) /n s . We note that for this relation to hold the time interval between samples St s needs to be 
somewhat larger than tq ~ 0.06, in order to ensure a number n s of statistically independent evaluations of j p , X y 
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Fluid 


Vpc 


T 


V 


Numerical 

((jp.xy) ) 


Theoretical 

(j\jp,xy) ) 


WCA 


81 


1 


1.75 


0.51 


0.60 


WCA 


173 


1 


1.75 


0.40 


0.41 


WCA 


138 


1 


1.75 


0.33 


0.38 


WCA 


338 


1 


1.75 


0.25 


0.24 


WCA 


2778 


1 


1.75 


0.08 


0.06 


WCA 


2778 


1 


1.75 


0.04 


0.03 


LJ 


121.5 


4.0 


2.12 


1.09 


1.08 (1.19) 


LJ 


121.5 


2.0 


1.90 


0.66 


0.72 (0.71) 


LJ 


121.5 


1.0 


1.75 


0.39 


0.49 (0.43) 



TABLE II: The variance of the microscopic momentum flux for p = 0.8. The theoretical results come from Eq. O using 
tg = 0.06; those in parentheses are derived from Eq. H16I . 



For the shear flows considered here the mean momentum flux can be expressed as 777 and the signal-to-noise ratio 
is E = 777/5. We now demand E > 0(1) to obtain conditions on Vpc and At av ; using Eq. i|T%|) one obtains 

E * = tf- T -E-^ av V PC > 0(1) (20) 
1 dt s 

and using 5t s > tq one gets, 

Vpc Atav > (21) 
TV 

Equations (|20|l and (|21|l give some control over the amplitude of the fluctuations of the coarse-grained (space-time 
averaged) stress tensor. For instance, in the case of steady state flows the signal-to-noise ratio can be increased by 
enlarging Vpc, Atav, or both. But if the flow is not homogeneous in space or if it is time-dependent, Vpc and At av 
are bounded above by the requirement of spatial and temporal flow resolution. In other words, to resolve a flow whose 
shortest characteristic time is w~^ x and whose maximum wavenumber is & max , the temporal and spatial resolution 
conditions are Afj „w~ ax < O(0. 1) and Axpck max < O(0. 1), where Vpc = ApcAxpc (see Fig. P). Inserting these 
conditions into Eq. (|21|l one gets that the range of frequencies and wavenumber that can be resolved is bounded 
above by ku; < ^rjT' 1 A~p 2 c . 



VI. APPLICATION TO OSCILLATORY SHEAR FLOW 



In order to test the applicability of the full hybrid scheme under unsteady flows, we consider the flow of an 
incompressible and isothermal fluid between two parallel walls in relative oscillatory motion. In our test flow, the 
simulation domain is < x < L x and it is periodic along y and z directions. The particle domain occupies the region 
x < lp, and it includes the LJ liquid and the atomistic wall composed of two layers LJ particles at x < 0. The 
continuum domain comprises the region x £ [lc,L x ]. The sizes of the simulation domains are within the nanoscale 
L x ~ 50a, and lp ~ 15a, while the width of the handshaking region are lp — lc ~ 5er. The flow is uniquely driven 
by the oscillatory motion of the x = L x wall along the y direction, meaning that the mean pressure is constant 
throughout the domain and there are no transfers of mean energy or mass in the x direction (perpendicular to the 
P^C surface). The mean flow carries transversal momentum by diffusion only, and the equation of motion for the 
y-velocity is du/dt — vd 2 u/ 'dx 2 , with boundary conditions u(0,t) = and u(L,t) = u wa ii(i) = M ma x sin(o;t). This 
equation can be solved analytically [a, 1241 |25j . The flow profile has a maximum amplitude at the moving wall and the 
momentum introduced by its motion penetrates into a fluid layer of width 8 ~ ^pnvjj . Beyond this layer (usually 
called the viscous layer) the flow amplitude tends to zero diffusively as it approaches the other wall held at rest. 
Therefore, the maximum shear rate attained inside or near the viscous layer is of order 7 ~ « max /i5. Inserting this 
relation into the signal-to- noise condition (|21|l . we find 

pu 2 max At av > nf- 1 , (22) 
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which clearly means that in order to attain a signal-to-noise ratio E larger than one, the mean kinetic energy of 
the flow integrated over the averaging time At av needs to be larger than the net energy due to fluctuations over a 
period of the mean flow. We note that at low enough frequencies, / > v/L 2 , diffusion is the fastest process in the 
system (momentum has enough time to spread mass over the whole domain) and in such cases the correct condition 
to guarantee E < 1 is given by Eq. (|21|1 . with 7 ~ w max /L 2 . 

In order to resolve the temporal variation of the flow the averaging time has to be a fraction of the shorter 
characteristic time of the flow, which corresponds here to the period of the oscillation, This condition, At av f > 
O(O.l), can be inserted into Eq. 11' I'D to obtain an estimation of the range of wall velocities that the scheme can 

1 /2 

resolve: u max > 5 ( pVpa ) ' ^ or = ^ = ^ an< ^ ^ PC = ®0-®®) tne above inequality yields u ma x > 0.5. 
Provided that At av f is fixed to a low value, we note that the former condition is independent of the external frequency 
because the shear stress is approximated by 7 ~ u max /5. This is only valid up to the end of the viscous layer because 
7 = 7(2;) becomes very small beyond this layer (x > 5). Also, using the analytical flow solution |g] it can be easily 
shown that for a fixed u max , the local shear rate within the viscous layer decreases slightly with the frequency. In 
fact, the characteristic Reynolds number of the flow goes like J -1 / 2 : Re* — u m3X S/u ~ w max / -1 / 2 ^ -1 / 2 . As an 
aside, the flow considered here (the oscillatory shear flow in between infinite planes or stream- wise periodic boundary 
conditions) is linearly stable |16|. As shown in a recent paper, flow instabilities in the form of three-dimensional rolls 
do arise at Re = u max L x /i> ~ 500 only if the zero mass-flux condition along the stream-wise direction is imposed by 
confinement of lateral walls ^(|. The range of Reynolds number considered here was such that Re < 200. 

In order to test the hybrid model and the above observations on the amplitude of fluctuations, we performed 
oscillatory shear simulations for values of u max above, near and below the threshold given by Eq. (|22Jl . In the 
following calculations the domain geometry is L x = 30 and in the periodic directions L y — L z were set to values 
between 6 and 9. Calculations for the large amplitude flow are in excellent agreement with the analytical solution as 
illustrated in Fig. where u max = 10, / = 0.01, and At av = 1. 

Figure shows several snapshots of the velocity profiles in the case of a large amplitude flow (u max = 10, / = 0.01) 
of a LJ fluid at density 0.8 and temperature T = 1.0. The P region spreads from the atomic wall at x = to x = 15 
and the C region comprises x = [11.8, 30]. The centres of the P— >C and C^P cells were set 3.2 apart and their width 
was 1.6. The outcome of the hybrid simulation is in very good agreement with the analytical solution, plotted in 
dashed lines. 

The calculations shown in Fig. [|)]were performed for a velocity coupling of a — 0.5. Figure compares this result 
with another calculation using a — 0. In the latter case, discontinuities of velocity are clearly observed within the 
overlapping region, its maximum size, around 0.35, being consistent with the size of the velocity (and momentum 
flux) fluctuations within the overlapping region. It is interesting to note that such discontinuities can be smoothed 
out by using values of a as small as 0.2. In other words, the model is not very sensitive to the actual value of a, 
provided a > 0. 

Some simulations with the same fluid (LJ at p = 0.8 and T — 1) but using a smaller flow amplitude, M max = 0.5 
are shown in Fig. [S] The averaging time was chosen to be At av = 10. The case of Fig. corresponds to a frequency 
of / = 0.01 and it lies approximately at the threshold given by Eq. I|22|l . The instantaneous velocity at the P->C cell 
indicates that the noise amplitude is nearly equal to the flow amplitude; as a consequence, the corresponding time- 
averaged velocity (shown in Fig. |SJj) wanders around the analytical curve, clearly showing the traces of instantaneous 
fluctuations. At smaller frequencies the momentum spreads over the whole domain and the maximum velocity attained 
in the neighbourhood of the overlapping region is larger (recall that the Reynolds number depends on Z" 1 ^ 2 ). This 
explains why the simulation with / = 0.002 in Fig. [Bp is in better agreement with the analytical trend than that 
for / = 0.01. Finally, in Fig. [SJ; we present a case with M max = 0.5 but at a larger temperature T = 4. The flow 
amplitude is then below the threshold given by Eq. 122(1 and forces arising from thermal fluctuations are larger than 
those arising from the mean flow. 

VII. CONCLUDING REMARKS 

We have presented a hybrid finite volume formulation for the P— >C coupling in our atomistic-continuum scheme 
based on the exchange of fluxes between both domains. The scheme has been used to perform hybrid simulations of 
the oscillatory shear flow driven by the harmonic motion of one of the walls of a nano channel. The particle domain (P) 
comprises one of the walls and a part of the adjacent fluid layer, while the remaining fluid within the slot is described 
by continuum fluid dynamics (CFD). The results show that the flux-coupling scheme is able to solve unsteady flows 
up to high characteristic frequencies corresponding to a Stokes number of St = 2irfL x /v < 0(100). 

A central part of this study has been focused on the limitations of the scheme due to the effect of the fluctuations 
arising from the particle domain. We found that in flows with low shear rates 7 < O(10~ 2 ) fluctuations of velocity 
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and momentum flux may affect the continuity at the overlapping region. For 7 < O(10 2 ) these fluctuations induce 
differences of the P and C velocities which are of the same size as the maximum flow velocity. This affects the 
flux-scheme because the momentum flux injected into the P region is measured from the local velocity gradient in the 
C domain. We were able to largely avoid this problem by implementing a finite volume scheme that uses a hybrid 
definition of the velocity gradient, constructed from a linear combination of the C- velocity and the averaged molecular 
velocity at the flux boundary condition for the C domain. The algorithm derived for time integration of the C velocity 
has an extra term that produces a seamless velocity profile by gently driving the C-velocity towards the averaged 
P-velocity. The strength of the coupling term is proportional to a coupling parameter a G (0, 1] and the scheme 
is not very sensitive to the value of a, provided that a > (we used a ~ [0.1 — 0.5]). It is important to mention 
that the velocity coupling term only acts on the C domain and has no direct effect on the particle dynamics. Such 
a procedure avoids the use of Maxwell daemons within the overlapping region, which appear in schemes that rely 
on variable-coupling (2||[27|]. Another attractive feature of the P^C coupling presented here is that it guarantees 
the velocity continuity while still preserving the balance of fluxes in the coarse-grained dynamics. Indeed, the use 
of hybrid gradients may be extended to any desired variable. For instance, in the case of diffusion of mass density, 
Flekk0y et al. 8] used a purely hybrid gradient (i.e., with a = 1 in the present notation) to couple random walkers 
with a finite difference description of Fickian diffusion. 

If the signal-to-noise ratio E of the momentum flux is smaller than one the flow cannot be described using a 
deterministic scheme for the continuum flow. We have analysed this flow field resolution limit by calculating E using 
a route based on the Green-Kubo formalism. Resolution conditions are written in terms of the shear rate in Eq. 
(|21[) and in terms of the oscillatory- wall velocity in Eq. (|22[) . Simulations performed near and below the derived 
resolution limit were shown to be consistent with the theoretical considerations. If the flow amplitude is smaller than 
the threshold given by Eqs. i|21|) or IL'l'l) the effect of fluctuations should also be taken into account by describing the 
C domain through a mesoscopic fluid model that solves the stochastic Navier-Stokes equations [28l |29| . In the case 
of Fickian diffusion, Alexander and coworkers used a hybrid particle-continuum algorithm to study the fluctuating 
hydrodynamic limit • They considered a hybrid coupling scheme with both deterministic and stochastic continuum 
solvers and found that, if the correct stochastic representation is used within C, the mean and variance of the density 
are correctly coupled within the overlapping region. The authors showed that a deterministic hydrodynamic solver 
guarantees the coupling of the mean density fields, although in this case the variance falls exponentially in the direction 
of the continuum domain. We observed similar results in our simulations. The generalisation of our hybrid scheme to 
largely fluctuating flows is feasible and the implementation of a stochastic continuum solver would not substantially 
alter the scheme presented here. 
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FIG. 1: The domain decomposition of the hybrid scheme, (a) displays the set-up used in the present work (periodic in y and z 
direction) indicating the particle domain (P), the continuum domain (C) and the handshaking domain comprised of the C^P 
and P-^C cells, whose surface area is A. The P domain is contained in x < xcp and includes the wall around x ~ formed by 
two layers of LJ atoms in a hexagonal lattice (black circles). The C domain is xpc < x < L x and includes a boundary condition 
at x = L x , which imposes u y (L x ,i) = u wa u(t) (moving wall). Open circles represent particles whose dynamics are uniquely 
determined by their corresponding intermolecular forces. Particles entering the C— >P cell (hatched circles) feel extra forces 
coming from the continuum domain, (b) shows three consecutive cells of the finite volume discretization of the C domain: H, 
W = H — 1 and E = H + 1. The index of the cell centres runs from H = 1 to H = M. For the first cell H = 1, wherease the 
west cell W(— 0) is outside the C domain and within the P— >C buffer, as indicated in (b) (see Sec. 11111 . The centre-to-centre 
distance is Ax e — xe ~ xh and the face-to-face distance is Axh — x e — x w . 
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FIG. 2: The y-velocity versus the x coordinate for tests done in a Couette flow with wall velocity M max = 0.5. The particle 
region comprises x < 20 and the continuum domain is x G [16.36,30]. The mean density is p = 0.8 and the temperature is 
T — 1.0. Vertical lines indicate the average x position of particle wall (x = 1), the P— >C interface xpc = 16.31 and the extent 
of the P— >C and C— >P cells, Axpc = 1.81<r. The upper figure shows the y-velocities averaged over ~ 800r obtained for a — 1 
and a = 0.2. The lower figure shows the y- velocities obtained using a = 0, averaged over ~ 700r. 
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FIG. 3: The effect of the a parameter in Eq. il l U on the continuity of the y-velocity at x\ = xpc + Ax/2 (xpc = 16.31 and 
Ax = 0.9) for the same Couette now tests shown in Fig. [2] From top to bottom the figures correspond to a = 0, a — 0.2 and 
a — 1. The continuum velocity at a; = xi is labelled by u y and the mean particle velocity within the region \x — xi| < 0.9 
is denoted as (v y ). Dotted lines correspond to velocities averaged over a short time interval At a v = It, and solid lines are 
the same velocities averaged over 20r. The velocities u y and (v y ) are indicated in (a), while in (b) and (c) their values are 
almost identical. Note that the size of the velocity discontinuity observed for a = is similar to the amplitude of the velocity 
fluctuations (dotted lines). 
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FIG. 4: The particle density versus the direction normal to the atomistic wall at x < for a WCA-LJ fluid of density p c — 0.8 
and temperature T = 1.0. The density profile obtained with the hybrid scheme is compared with the outcome of a full MD 
simulation with a second atomistic wall at x — 50a. (a) corresponds to po = 0.7 in Eq. H14|l and (b) to po = 0.5. The width 
of the C^P cell is Axcp — 2.3a. The locations of the P— >C and C^P cells are also shown (see Sec. If VB I for further details). 
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FIG. 5: The velocity within the overlapping region for an oscillatory flow with wall velocity u ma ii(t) = Umax sin (27r/t), of 
amplitude tt max = 10 and frequency / = 0.01. The particle domain comprises x = [0,21] and the continuum x £ [17,30]; the 
extent of the periodic directions is L y — L z — 9. The nanochannel is filled with a LJ fluid at p — 0.8 and T = 1.0. We plot the 
instantaneous particle velocity at the P—+C cell (noisy signal) and the time-averaged velocity over At av = 1 at the C^P cell. 
Dashed lines correspond to the analytical solution of the flow. 
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FIG. 6: Snapshots of the j/-velocity versus the x coordinate at several instants in time for the same flow as in Fig. The 
particle domain comprises x € [0,15] and the continuum region x 6 [12,30]. The mean velocities within the P region are 
calculated within slices of width 1.6 and averaged over At av — 1; these velocities are plotted using different symbols for each 
time instant. We used a = 0.5 in Eq. 1111 . The solid lines correspond to the solution within the continuum region, solved via 
a finite volume scheme, and the dashed lines are the analytical solutions. All quantities are given in reduced LJ units. 
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FIG. 7: Comparison of the velocity profiles obtained with a = 0.5 and a = in Eq. il It . The parameters are those of Figs. 
|S|and|S| Results were obtained with At av = 1. The velocity discontinuity observed in the a — calculation (about 0.4) is of 
the same order of magnitude as the variance of the mean velocity within the P— >C cell. All quantities are given in reduced LJ 
units. 
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FIG. 8: The mean particle velocity in oscillatory shear flows with tt m ax = 0.5 and p — 0.8. The molecular dynamics domain 
comprises x € [0, 20] and the continuum x € [16.36, 30], the velocity is evaluated within the P^C cell (around x = 17.27). (a) 
corresponds to T = 1 and / = 0.01; (b) to T = 1 and / = 0.002 and (c) to T = 4 and / = 0.01. Solid lines correspond to 
the time-averaged mean particle velocity (in time intervals At a v = 10), dotted lines to the instantaneous mean velocity and 
dashed lines to the analytical solution of the flow. In (a) we show the time-averaged mean particle velocity {v)p along with 
the continuum velocity u at the same location. 



